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Abstract. We present a finite difference version of the eth formalism, which 
allows use of tensor fields in spherical coordinates in a manner which avoids polar 

~ singularities. The method employs two overlapping stereographic coordinate patches, 

with interpolations between the patches in the regions of overlap. It provides a 

new and effective computational tool for dealing with a wide variety of systems in 

<~j which spherical coordinates are natural, such as the generation of radiation from an 

isolated source. We test the formalism with the evolution of waves in three spatial 
dimensions and the calculation of the curvature scalar of arbitrarily curved geometries 

(T) on topologically spherical manifolds. The formalism is applied to the solution of the 

Robinson- Trautman equation and reveals some new features of gravitational waveforms 

CN in the nonlinear regime. 
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1. Introduction 

The inspiral and merger of a binary black hole system is anticipated to be the prime 
source of gravitational waves for future wave detectors. Calculation of the emitted 
waveform, by means of emerging supercomputer technology, is the goal of the Binary 
Black Hole Grand Challenge [1]. In this paper we present a formalism which plays a 
strategic part in extracting the waveform produced in this three dimensional problem 

[2]- 

At large distances from a compact source the wavefronts of any radiation held 

become spherical; this leads naturally to the use of a spherical coordinate system. 

Indeed, spherical coordinates and spherical harmonics are standard analytic tools in 

the description of radiation. But the use of spherical coordinates in numerical work 

leads invariably to the vexing problem of coordinate singularities at the origin and 

along the polar axis. Finite difference approximations suffer particularly because they 

have no natural way of enforcing the correct boundary conditions on their solutions. 

As a result, spherical coordinates have mainly been used in axisymmetric systems, 
where the polar singularities may be regularized by standard tricks. In the absence 
of symmetry, these techniques do not easily generalize and they would be especially 
prohibitive to develop for tensor fields. 

This paper presents an approach that should be of interest to researchers working 
in numerical relativity in problems where spherical coordinates are a natural tool. 
The eth formalism [3, 4] and the associated spin-weighted spherical harmonics [4, 5] 
allow a simple and unified description of vector and tensor fields, without the undue 
complexity of traditional vector and tensor harmonics. Unfortunately, these techniques 
have remained obscure and unfamiliar to many researchers who would benefit from 
them. 

We present here a finite difference version of the eth formalism, with applications 
to numerical relativity. The framework for the computational application of the eth 
formalism is developed in Sec. 2. Although our presentation is on the algorithmic 
level, we present enough detail to direct a well defined transition from algorithm to 
code development. In Sec. 3, we present tests and applications of the formalism which 
demonstrate both its accuracy and usefulness. The hrst application, the evolution of 
scalar waves is illustrative of the techniques necessary to carry out a 3-dimensional 
evolution using spherical coordinates. The second application, the numerical calculation 
of the curvature scalar of a topologically spherical manifold illustrates the global and 
tensorial features of the method. This calculation is at the heart of several important 
problems, such as the computation of the Hawking mass [6]. The third application, 
solution of the Robinson- Trautman equation, reveals new and unexpected nonlinear 
properties of the gravitational waveform. This work illustrates how computational 



relativity can greatly benefit from approaches based upon the best analytic tools. 

2. The Eth Formalism 

We present here a finite difference treatment of the sphere, based upon (i) the standard 
method of describing the global differentiability of functions by means of local coordinate 
patches, and (ii)the eth formalism for describing differentiation of fields on the sphere. 

2.1. Patching 

Let 6 and </> label the points on the sphere. Then, the real and imaginary parts of 
the stereographic coordinate (n = tan(#/2)e 8?i provide a smooth coordinatization of the 
sphere excluding the point 6 = ir at the south pole. Similarly (s = 1/CiV provides a 
smooth coordinatization except at the north pole. These two stereographic coordinate 
patches are sufficient to cover the sphere. It would be possible, in principle, to cover 
the sphere with a number of non-singular (#, </>) patches, by selecting several directions 
(6i } (f>i) on the sphere and associating with them a local coordinate patch, and then 
excluding from those patches the (local) polar region. Stereographic coordinates, provide 
the most economical choice, where only two non-singular patches are used. 

A smooth scalar held *& on the sphere may be represented in terms of smooth 
functions ^ s of (s on the lower hemisphere \(s\ < 1 and as smooth functions ^jv of 
(n on the upper hemisphere |(jv| < 1- The continuity of *& is ensured by requiring 
^slCs) = ^n(Cn) at the equator but in order to ensure smoothness this equality must 
extend to a common overlap region about the equator. Our hrst task is to implement 
this in a finite difference treatment. 

In order to construct the grid, let (S 1 , S 2 ) and (N 1 , N 2 ) be the real and imaginary 
parts of (s and (jv, respectively. Let the grid consist of the points S l = Si/M and 
N l = rii/M, where s 8 - and n 8 - are each pairs of integers with —M<Si<M and 
— M<rii< M. Thus the grid points at the intersection of the real and imaginary 
axes with the equator correspond to s 8 - = (M, 0) and s 8 - = (0,M); and n 8 - = (M, 0) 
and rii = (0,M). The grid points at the poles correspond to s 8 - = and n 8 - = 0. Note 
that our choice of rectangular grid domains is not unique [7]. Other strategies can be 
followed, to economize the number of points in the overlap region between the S and 
iV grids, but the rectangular choice lends itself to a simple implementation at the grid 
boundaries. 

Given a smooth function ^ on the sphere, we now give instructions for representing 
it on the grid. Denote the values of V I / 5(C5) at points on the S'-grid by \&s SliS2 ; and the 
values of ^n(Cn) at points on the iV-grid by x ^Nn 1 ,n 2 - I n the overlap region, the grid 
representation must be consistent with the condition V I / 5(C5) = ^n(1/Cs)- However, 
only in exceptional cases will a point that lies on the S'-grid also lie on the iV-grid. To 



deal with this, we introduce 

_ (rai -m 2 )M 

which gives the value of (s at a point which lies on the iV-grid. Similarly, we define 

(s 1 -is 2 )M 
Cns Si , S2 - {s2 + 4) • (2) 

In the same way, in the overlap region, ^5 and ^jv determine the values x $ss 1 ,s 2 and 
x &NSs 1 ,s 2 at points which lie on the S'-grid; and ^Nn u n 2 and ^ SNn u n 2 at points on the N- 
grid. However, the grid values x $ss 1 ,s 2 are weaker information than the smooth function 
^ s and do not by themselves determine the value of *& at points on the iV-grid. Similarly, 
the values ^Nn u n 2 do not by themselves determine ^ at a point on the S'-grid. In these 
cases interpolation is necessary to evaluate functions on one grid at points lying on the 
other grid. 

In order to see how this applies to differentiation, consider the gradient of ^ , with 
components ds'^s(Cs) in the S'-patch; and <9jv° v I / jv(Civ) i n the iV-patch. Then, at points 
which do not lie on the grid boundary, the partial derivatives can be approximated by 
centered finite differences in the standard way. For example, at the point S l = s 8 A, the 
appropriate finite difference approximation is 

d sl q> s = ^£fi±Mi_ l^izhfL + (A 2 ). (3) 

For a point on the grid boundary, e.g. the point S 1 = (I,s 2 A), we use the same 
approximation except that the held at the virtual grid point s 8 - = (M + I,s 2 ) is 
approximated by the value x $ns(m+i),s 2 obtained from interpolating ^n- In order to 
achieve an error of second order in the derivative the interpolation error must be of 
third order. 

This provides second order accurate finite difference representations of the gradient 
of ^ at all grid points and, in fact, two representations, corresponding to ds^s(Cs) 
and <9jv° v I / jv (Civ)? at points which lie in the overlap between the two grids. The analytic 
expressions are related by the transformation 

d< s ^s(Cs) = ^d (N V N (C N ). (4) 

cJ Cs 

As a check on the finite difference approximation, the finite difference equivalent of (4) 
must hold to second order accuracy. Equivalently, the gradient may be converted into a 
scalar held by taking basis components and the scalar held compared between patches. 
This is the strategy of the eth formalism that we shall consider later. 

Second derivatives of a scalar held with respect to the coordinates may also be 
approximated to second order by central differencing. The interior points in each grid 



patch are handled by the standard techniques and the points on the grid boundary 
are handled by introducing a virtual grid point, as above. The following fourth order 
interpolation scheme ensures that both the hrst and second derivatives are approximated 
to second order accuracy at boundary points. 

First, we extend the grid by one point on each coordinate direction. The points at 
the S grid with coordinates s 8 - = (M + l,s2) will define the computational rightmost 
boundary of the grid, etc.. Centered derivatives in the S\ direction require the held 
value at the virtual grid point s 8 - = (M + 2, s2). Note that this virtual point maps into 
a point on the N grid which has at least two neighbors on each direction who also fall 
inside the computational boundary of the N grid, namely their grid coordinates satisfy 
the relations |n 8 | < \M + 1|. 

We interpolate the value at the virtual point by the following procedure. First, we 
evaluate the coordinates n 8 - of the virtual point, which according to (2) are given by 
rii = S\M I (s\ + s^), n 2 = —S2M I (s\ + S2). We then compute the indices (ib,jb) of 
the lower left corner of the grid cell on the N grid where the virtual point falls, i.e. 
ibA < rii < (ib + 1)A and j&A < n 2 < (jb + 1)A. The nine grid cells whose corners 
are given by the points (i,j) with ib — I < i < ib + 2 and jb — 1 < j ' < jb + 2 fall 
entirely inside the N grid. We then fit a cubic Lagrange polynomial to each of the four 
set of points obtained by keeping the j index fixed, and obtain the value of the held 
at the points (si, j), J = jb — 1, ..., jb + 2. Another cubic Lagrange polynomial is then 
fit through these four points to yield the value at the virtual point (si,s 2 ). Note that 
all grid points at the boundary can be treated with the same algorithm. This method 
has the computational advantage that it vectorizes on machines with the appropriate 
hardware, such as Cray supercomputers. 

We have carried out checks of the interpolation error. Even for very coarse grids, 
the measured convergence rate is ss 3.9, in good agreement with the theoretical limit of 
four. 

2.2. Eth 

The eth (3) formalism gives a compact and efficient manner of treating vector and tensor 
fields on the sphere, as well as their covariant derivatives. In addition, the associated 
spin-weighted spherical harmonics provide a very tidy alternative to vector and tensor 
spherical harmonics. We give here a brief description of the basic ideas and the details 
of how it works in some simple cases. 

The line element for the distance between neighboring points on the sphere is given 
by the quadratic form 

ds = q a bdx a dx (5) 

where the components q a b of the metric tensor depend upon the choice of coordinates. 
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In standard spherical coordinates this takes the form ds 2 = dO 2 + sin 2 Odcj) 2 but here we 
use stereographic coordinates for which 

ds 2 = 4(1 + CsCsr'KdS 1 ) 2 + (dS 2 ) 2 ] (6) 

in the North patch, with a similar expression holding in the South patch. The eth 
formalism represents the metric in terms of a complex basis vector q a , 

q a b = (q a qb + q a qb)/2, (7) 

where q a q a = and q a q a = 2. (Note that this departs from other conventions [4] to 
avoid unnecessary factors of \/2 which would be awkward in numerical work.) Since the 
real and imaginary parts of q a are both vector fields of unit length, this basis cannot be 
assigned in a smooth manner over the entire sphere so that two patches are necessary. In 
the S patch we make the choice q§ = (l + £s£s)((^ + ^2)/2, so that its real and imaginary 
parts line up with the S axes. Similarly, in the N patch, q^ = (1 + CnCn)($i + ^2)/^- 
Any two complex bases q a and q a are related by a unitary transformation q a = e ta q a} 
where the phase a is a real valued function. In particular, in the overlap between the 
patches, q^ = e lCi q a s , with e lCi = —(s/(s- 

Introduction of this complex basis allows us to represent any vector held on the 
sphere, with components v a , in terms of the complex scalar held v = q a v a . Here it is 
implicit that v represents either vs = qs v a or v^ = q%v a , depending upon which patch 
is being used, with v^ = e ta vs in the overlap region. A held v with this transformation 
property is called a spin-weight 1 held. 

Similarly, a tensor held on the sphere, with components w a f, } can be represented by 
scalar fields. First we decompose w a f, into its symmetric-trace-free part, its trace part 
and its antisymmetric part: 

w a b = t ab + -q ab + i-(q a q b -q a q b ), (8) 

where p = q cd w C( i and u = i(q c q d — q c q d )w c d/2. The (real) scalar fields p and u are 
independent of choice of basis and are called spin-weight zero fields. The symmetric, 
trace free tensor held t a b can be represented by the complex scalar held t = t a \ ) q a q h . 
Under change of basis we then have tjy = e 2ta ts and t is called a spin- weight 2 held. Thus 
an arbitrary tensor held can be represented by two real spin- weight fields and a complex 
spin-weight 2 held. These spin-weighted fields are the irreducible representations (of 
the unitary group of basis transformations) contained in the tensor held t ab . Note that 
i = tabq a q h alternatively represents the trace-free symmetric part of w ab as a held with 
spin- weight minus 2. 

2.3. Derivatives of First and Second Order 

In treating derivatives of tensor fields, it is again convenient to represent them in terms 
of spin-weighted scalar fields. This can be accomplished by taking basis components of 
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covariant derivatives V a . In the case of a scalar field, V„f = d a ^ , which we represent 
by the spin- weight one quantity 3^ = q a d a x $ } which represents the components of the 
gradient in terms of a complex held. Thus, at the grid point s 8 - in the S patch, centered 
approximations to the derivatives give 

Ss* 5 = i±lf|±iI^!($ Ssi+1)S2 _ tf Ssi _ M2 + t^ SsuS2+1 - ^ Ssi , S2 -i) + 0(A 2 ) (9) 

The overlap condition between patches is Qj^^n = e ta (3s^ s- 

In the case of a vector held, the covariant derivative, with respect to the unit sphere 
metric, is 

V a v b = 8 a v b - F ab v c (10) 

where the connection coefficients are determined from the metric by 

T^ = q cd {d a q db + 8 b q ad - d d q ab )/2. (11) 

This is equivalent to the spin-weight held Bv = q a q b V a v b and the spin-weight 2 held 
<5v = q a q b V a v b . Here, as above, we have v = q a v a . These complex fields contain all the 
information about the gradient of v a . Finite difference approximations can be obtained 
by hrst reexpressing 

Qv = q a d a v - Tv (12) 



and 



where 



Qv = q a d a v + Tv, (13) 



T = -q a q b V a q b /2 = (. (14) 

For a second order accurate finite difference approximation, the hrst term in either (12) 
or (13) may be approximated by a centered difference, as in (9), and the second term 
evaluated using the grid values of ( in (14). This leads to finite difference expressions 
in both the S and N patches, satisfying 3jvUjv = ^s^s and 3jvUjv = e 2ta &svs in the 
overlap. 

To treat differentiation of two-index tensor fields, it suffices to consider trace-free 
symmetric tensors represented by the pure spin-weight 2 held t = t ab q a q b . Then 
the covariant derivative of t ab can be expressed in terms of the spin-weight 1 held 
5t = q a q b q c, Vatbc and the spin-weight 3 held 5t = q a q b q c X7 a t bc . These can be reexpressed 
as 

8t = q a d a t - 2ft (15) 

and 

at = q a dj + 2rt, , (16) 



8 

with T given by (14). As before, these expressions have straightforward finite difference 
approximations in each patch, with the overlap relationships 3jytjv = e ta (3sts and 
5jvtjv = e 3ta (3sts- This procedure generalizes to treat derivatives of tensor fields with 
an arbitrary number of indices. 

Second covariant derivatives can be treated as products of first derivatives. However, 
this can lead to inefficient finite difference approximations. For example, the Laplacian 
may be written as 

D 2 ty = q ab V a V b ty = q a q b V a V b ty = 55*. (17) 

A straightforward finite difference approximation to 3 followed by an approximation to 
3 would involve next nearest neighbors, whereas a central difference approximation to 
the second derivative would involve only nearest neighbors. Both approximations are 
second order accurate but the latter is clearly preferable. There are various other second 
derivative expressions in which the latter approach should also be invoked. 

One such expression, which commonly occurs in dealing with a spin-weight 2 held t, 
is the spin- weight zero held 3 2 t. We will use this as an example of how such expressions 
can be approximated to second order using only nearest neighbors. In our choice of 
stereographic coordinates, we have 

B 2 t = q a q b d a d b t-2((m + (t). (18) 

Thus calculation of 3 2 t reduces to the previous calculation of <3t and the mixed 
derivatives d a d b t } which can all be approximated to second order by centered differences 
involving only nearest neighbors. At the grid boundary, this involves values at virtual 
grid points which can be determined to fourth order accuracy by the interpolation 
scheme. 

Two general operators expressing arbitrary hrst or second order derivatives acting 
on any spin-weighted held are defined here. Those operators have as arguments the 
held values over the sphere, the spin of the held and a set of binary values that indicate 
the nature of the operator. These operators allow for easy and versatile numerical 
implementation and cover the vast majority of expressions we expect to encounter in 
numerical relativistic calculations. 

The hrst order operator is defined as 

Di(s, e)V> = -(1 + x 2 + y 2 )(4V + ied y il>) + s(ex + iy)^, (19) 

where ( = x + iy, s is the spin- weight of ip and e = 1 gives the 3 operator while e = — 1 
gives the 3 operator. 

Similarly the second order operator is defined by 

D 2 (s, ei, e 2 )^ = {1(1 + x 2 + y 2 f[d 2 x - e 1 e 2 d 2 + t{e 1 + e 2 )d x d y ] 



+ -(1 + x 2 + y 2 )[((l + e l£2 + s(ei + e 2 ))x + i{2s + t x + e 2 )y)& 

+ ((2s + ei + e 2 )x + «(1 + e 1 e 2 + s(ei + e 2 ))y)ieie 2 d 2/ ] 

+ s [^( £ 2 - £ i) + {seie 2 + -(ei + e 2 ))(x 2 - e 1 e 2 y 2 ) 

+ (l + € 1 € 2 + s(e 1 + € 2 ))ixy]}ip. (20) 



The different combinations of ei, e 2 generate the following second order operators: 

dm, i) = a 2 , dm, -l) = aa, d 2 (m, i) = aa, z> 2 (-i, -i) = a 2 . 



3. Applications and Tests 

3.1. Wave Evolution 

An important hrst application and test of the approximation scheme developed in 
the previous section is the numerical solution of the 3-D wave equation in spherical 
coordinates. We start with the the discretization of the Laplace operator on the sphere. 
In terms of the complex coordinate (, we have 

D 2 m = D 2 (0, 1, -1)* = (1 + (0 2 d ( d^. (21) 

Since D 2 is a scalar operator, (21) gives the same value in the overlap region when 
evaluated in either the S or N coordinates. 

In retarded time coordinates the wave equation □ $ = takes the form 

D 2 g . . 

2y,„r - g, rr H — = 0, (22) 

where g = r$. We solve (22) by an explicit marching algorithm, from r = to null 
infinity [8]. 

The convergence and stability of the marching algorithm is dictated by the 
Courant-Friedrichs-Lewy (CFL) condition (which requires that the numerical domain 
of dependence includes the physical domain of dependence). In the present case, this 
requirement is strongest at the vertex of the outgoing cones. It limits the timestep by 

Au < kArA Sl As 2 , (23) 

where A; is a number of order one, related to the details of the startup procedure near 
the vertex r = 0. To good accuracy, our numerical investigations give the value k = 1 
for stable evolution. 

A good test of the patching scheme focuses on the accuracy of the numerical 
evolution. The linearity of the problem provides a useful family of exact solutions 

G lm (u,r, C,C) = , , 1 ^ 1/ r+ ,\ , ^ Y lm ((,(). (24) 



+ l)'+ 1 (u + 2r + l)'+ ll1m(C ' C) - 
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Expressions for the standard and spin-weighted spherical harmonics in stereographic 
coordinates can be found in [9, 4]. The characteristic scheme evolves a global solution 
to the wave equation and thus it is possible to check the global numerical error. The l 2 
norm of the error, for the evolution of initial data given by an / = 2, m = 2 harmonic, 
converges to zero with a measured 2.02 convergence rate. 

Figure 1 shows a sequence of snapshots of the radiation patterns, as seen at null 
infinity, for initial data which consist of a localized pulse in the North patch traveling 
towards the origin. Notice the smooth propagation of the held across the coordinate 
patches as the evolution progresses. 

3.2. Calculation of Scalar Curvature 

As a tensorial illustration of the forgoing methods, we consider a problem which arises 
in many different contexts in general relativity: Given the metric h ab of a topological 
sphere, calculate the scalar curvature. There are many ways to formulate this problem 
analytically. Here we present a treatment, in terms of an auxiliary unit sphere metric, 
which is flexible enough to be of broad usefulness in numerical applications. 

The metric is uniquely determined by its unit sphere dyad components K = 
h a b<l a <t/2 an d J = h ab q a q b /2, which are of spin weight zero and two respectively. The 
dyad components of the inverse metric are h ab q a q b = 2K/H and h ab q a q b = —2J/H, where 
H = (K 2 — J J) = det(h ab ) / det(q ab ) . Also note that in our conventions the alternating 
tensor of the unit 2-sphere is e ab = iq[ a q b ] and (33 — 55)?y = 2sr/ for a spin-weight s field 

V- 

Let D a represent the covariant derivative associated with h ab [10]. Then its curvature 

determines the commutator 

(D a D b - D b D a )q cd = h ef (R abcJ q ed + R abd jq ce ), (25) 

where, in two dimensions, R a bcd = Rh a [ c hd]b- We let the tensor field 

C c ab = \h cd (V a h db + V b h ad - V d h ab ) (26) 

represent the difference between the connection associated with D a and the unit sphere 
connection, e.g. (D a — \ 7 a )vb = ~C^ b v c . Then by contracting (25) with q a q b q c q d , we 
obtain an expression for the curvature scalar in terms of the dyad components of C c ahl 

RJ = 3C -QB + -(CA + AB-B 2 -BC), (27) 

where A = q a q b q c C c ab , B = q a q b q c C c ab and C = q a q b q c C c ab . 

Next, we use (26) to relate these to the metric components, 

A = ^(2KQK - KdJ - ,m,J) 
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b = ^(Kdj - ,m,J) 

C = ^-(KQJ-2JQK + JQJ). (28) 

Substitution into (27) then yields the expression (29) for the curvature scalar, after 
canceling a common factor of J from both sides: 

R = ^ 77 {2K - QQK + 5 2 J + ^-[2(QK)(KM - KQJ - JQJ) 

+ (aj)(jgj + Agj) + (gj)(jgj_ Agj)]} 

+ c.c. (29) 

where H = (K 2 — J J). The derivation of (29) is quite laborious and it would 
be extremely useful to develop symbolic techniques to carry out general dyad 
decompositions and reduce covariant derivatives to eth operations. In the form (29), the 
curvature scalar may be computed to second order accuracy using the finite difference 
techniques presented above. 

In the particular case where det(h a b) = det(q a b) } as in the Bondi [11] formalism, 
H = 1 and the metric is uniquely determined by J. In this choice of gauge the curvature 
formula simplifies considerably, 

R = 2K - mK + - (8 2 J + a 2 J) + 4r (3J3 J - 3-/3J) • (30) 

Note that in (30) the complex conjugation has been explicitly carried out. 

As a test, we choose J to be a simple spin-weight 2 harmonic, e.g., J = c(3 2 Y 2 o } which 
facilitates the calculation of exact expressions for the curvature. Scaling the metric data 
with the parameter c generates arbitrarily non linear curvature expressions which are 
used to check the numerical implementation of (30). 

An interesting check of the numerical procedure utilizes the global geometrical 
properties of curvature on spherical surfaces. The Gauss-Bonnet theorem for spherical 
topology implies 



1= <b RdS = Sir. (31) 

A global test of the computation of the scalar curvature can thus be performed by 
integrating the numerically evaluated curvature over the entire sphere. The surface 
area element is dS = 2iH 1 ' 2 (l + (()~ 2 d(d(. The integration must take into account 
the overlap between coordinate patches. A simple and natural choice is to use the 
equator (( = 1 as the smooth and symmetric boundary of the integration within each 
patch. A second order integration scheme in two dimensions is the simple four point 
center average method. We must pay careful attention to the boundary cells, i.e. those 
which overlap the equator. Their number scales like 1/A and if we were to omit their 



12 

contribution to the integral the order of accuracy would be reduced by one. To take 
them into account, the interior area of each boundary cell is approximated assuming 
that the equator looks like a straight grid line locally. The centered four point average 
of the cell is then weighted by its interior area divided by A 2 . 

Figure 2 presents the Gauss-Bonnet global test and confirms the convergence rate 
of 2. 

3.3. Robinson- Trautman Solutions 

Robinson- Trautman spacetimes [12] contain purely outgoing gravitational radiation, 
which decays to leave a Schwarzschild-like horizon [13, 14]. The Bondi news function 
for these spacetimes is given by [15] 

N = l - W -1 a 2 W, (32) 

where W(u, (*, C) satisfies the Robinson- Trautman equation 

i2d„w = w 3 (a 2 wa 2 w - wa 2 a 2 w). (33) 

The real and imaginary parts of the news function represent the two independent 
polarization modes of the radiation. The Bondi mass is given by the solid angle integral 
M = ( l/47r) § W~ 3 o?0. The Schwarzschild spacetime (for a unit mass black hole) 
corresponds to the solution W = 1. The linearized perturbation W = 1 + A(u)Yi m 
decays exponentially according to A = A(0)e~ u£ ( e+1 ^ £2+e ~ 2 ^ 12 . 

Accurate numerical evolution of a fourth order parabolic equation, such as (33), by 
means of an explicit finite difference scheme is a challenge because the C FL condition 
requires that the time step Au scale as the fourth power of the spatial grid size. 
Nevertheless, we have applied the 3 formalism to obtain an efficient, second-order 
accurate evolution algorithm, based upon a three time level Adams-Bashford [16] scheme 
with predictor (W) given by 

A?y 
W(u + Au) = W(u) + — d u [3W(u) - W(u - Au)], (34) 

and corrector 

Au 
W(u + Au) = W(u) + — d u [W(u) + W(u + Au)] + 0(Au 3 ). (35) 

Here the d u terms are to be considered reexpressed in terms of 3 operators using (33). 
The second order convergence of numerical solutions was confirmed in the perturbative 
regime using solutions of the linearized equation. 

This algorithm reveals some new and unanticipated features of the gravitational 
waveforms in the nonlinear regime. Figures 3 and 4 illustrate the Bondi news function 



13 

for the Robinson- Trautman space-time with initial data in the odd parity £ = m = 3 
mode 

W| u=0 = l + AK[y 3 3], (36) 

for the case A = 0.89. In order to provide a physically intrinsic time dependence, the 
news function is plotted versus Bondi time ug, which is related to the coordinate time 
u by dus/du = W. 

Figure 3 shows surface plots at representative times of the nonlinear contribution 
to the news function, i.e. AN = N(X) — d\N\\ =0 \. (The plots show the real part of 
AN in the south patch). The smoothness of the plots indicate the effectiveness of the 
algorithm in handling the high spherical harmonics generated by the nonlinearity. The 
highest harmonics damp quickly and the system asymptotically approaches a monopole- 
dipole combination at late times. 

Figure 4 displays the time evolution of the news function at a representative point 
on the equator. The graphs reveal an oscillatory behaviour qualitatively quite different 
from the pure exponential decay of the linearized solution. In addition, the behavior of 
the two polarization modes is quite different. Approximately 5% of the initial mass is 
radiated away during this simulation, which ends at the time of black hole formation. 
This is in the range of energy expected to be radiated during the inspiral of a binary 
black hole system. In this regime, we have found (by varying the value of A) that these 
oscillatory waveforms cannot be well approximated even by second order perturbation 
theory, which emphasizes the importance of an accurate computational treatment. 
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Figure 1. A sequence of radiation patterns at infinity for initial data consisting of 
a localized pulse which is traveling towards the origin. The left column represents 
the north patch, the right column the south patch. The evolution proceeds from 
top to bottom. The first pair is at retarded time u = 0.04 and displays the initial 
radiation resulting from scattering of the pulse off the centrifugal barrier created by 
its high angular momentum (/ multipole values). This radiation shows up in the same 
hemisphere as the initial data. At the next two retarded times, u = 0.16 and u = 0.32, 
radiation appears in all directions as the pulse crosses the origin defining retarded time. 
Finally, at u = 0.8 the pulse has completely passed by the origin and the radiation 
ceases. The grid size is M = 18, i.e. 37x37 points per patch. 
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Figure 2. The global Gauss-Bonnet convergence test. The figure displays the 
convergence of the numerically evaluated integral I = § Rds to the exact value of 
87r. The seed for the metric data is an / = 2, m = harmonic and it is scaled over 
six orders of magnitude. The convergence is uniformly second order. Note that for 
very low amplitudes the error is dominated by the surface integration procedure and 
is independent of the data. 
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Figure 3. Comparison of linear versus non-linear evolution of outgoing radiation 
fields. The graphs show the difference between the Bondi News for the perturbative 
and the numerically evolved spacetimes for (a) ub = 0., (b) ub = 0.07, (c) ub = 0.14 
and (d) ub = 0.25. At early times («g = 0.07), a complex angular structure reveals 
the presence of high spherical harmonics which damp in a short period of time. At 
later times («g = 0.25), the system approaches a monopole-dipole combination. 
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Figure 4. The Bondi News function at the point (6, <f>) = (ir/2,Tr/2), for initial 
data corresponding to Y33. Both the real (a) and imaginary part (b) are shown. A 
perturbative calculation to first order predicts that the real part at this point is zero, 
with the imaginary part decaying exponentially. In contrast, the real part of the 
non-perturbative News is markedly different from zero, and both components show 
oscillations at early times. 



